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ABSTRACT 



A stable numerical procedure for solving the coupled, 
nonlinear equations of electron and ion conservation and 
electrical potential (Poisson's equation), is described, A 
two-dimensional model with periodic active sites on a flat 
plate is utilized to obtain both qualitative and quantita- 
tive results which clearly illustrate the self-generating 
sheath and ambipolar regions adjacent to a non-emitting 
electrode. The active sites are characterized as voltage 
sources and by electron densities depressed from both the 
non-active wall and the free-stream values. Various cases 
of species density at the nonactive wall are evaluated. 
P.ecombination/ion izat ion plays an important role in estab- 
lishing the boundary layer behavior in the ambipolar region 
and in the dimensionality of the problem formulation. 
Application is for Nitrogen gas at one amagat in a pulsed 
discharge of about 50 micro-seconds in duration. 

Results clearly demonstrate the boundary layer nature of 
the species density at the electrode. A sheath length of 38 
microns about a 35-volt active site and 5-20 microns along 
the non-active portion of the electrode is established. 

Joule heating is determined to be not important. 
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I. INTRODUCTION 

A. BACKGROUND 

Losses at the electrode boundary regions play an impor- 
tant role in high energy density molecular gas flow lasers. 
Of particular interest to this study is the Carbon-Dioxide 
glow-discharge plasma, of which Nitrogen is the principal 
energy bearing constituent. 

The main energy loss mechanisms in a collision dominated 
plasma are found in the sheath and ambipolar regions. The 
sheath voltage loss occurs directly as a result of Debye 
shielding, which creates a space charge layer abutting an 
electrode or wall. The ambipolar losses are associated with 
ambipolar diffusion, relating to the charged species density 
being less than the free-stream or equilibrium density. 

This loss exists even in the absence of current flow. 

An initial premise of probe theories that takes account 
of collisions between charged particles, atoms and molec- 
ules, is that the plasma is quasi-neutral up to a certain 
distance from the probe and that the electric fields are 
localized within a narrow sheath region adjacent to the sur- 
face of the probe. The problem of making the transition 
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from the quasi-aeatral solation to the sheath solution was 
usually resolved through matching techniques. The descrip- 
tion of a collision dominated sheath together with the qua- 
si-neutral ambipolar region presents an imposing set of cou- 
pled nonlinear partial differential equations, for which 
there exists no known exact solution. This study reports a 
numerical method that describes the sheath and ambipolar 
regions in a self-consistent fashion. 

In the hydrodynamic approximation (i.e., collision domi- 
nated), the dimensions of the space charge region (the 

sheath length) is larger than the mean free path '• / " of the 
species particles. Additionally the sheath characteristic 
length is very much smaller than the characteristic length 
of the plasma boundary layer. This puts our study in the 
collisionally thin sheath region, •' Sef . 1]. 

However, the plasma is not considered to be in thermal 
equilibrium. In fact for plasmas in thermal nonequilibrium, 
the electron temperature is governed by the electron energy 
equation, which includes the effects of the electron energy- 
conduction and electric field as well as collisional energy 
exchange between electrons and heavy-gas particles. The 
consideration of negligible energy exchange during colli- 
sions (thermally frozen) does not generally result in a 
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constant electron temperature (T^) ; particularly in the 
sheath region, even if other properties such as the ion 
temperature (I^) are constant. Therefore the electron 
energy equation is necessary even in the case of the ther- 
mally frozen situation. 

An ionizer/sust ainer discharge pumped laser requires 
that relatively high E/n values be distributed as uniformly 
as possible throughout the bulk of the gas medium. The 
electrode regions however represent a substantially higher 
and more inhomogeneous B/n which may impose an instability 
on the remainder of the discharge. In the region of the 
electrode, the B/n is significantly higher than the undis- 
turbed plasma. In this region the elevated electron temper- 
ature, which is a direct function of the sustainer field, 
will increase ionization, which could lead to arcing and 
breakdown. 

In the so-called vacuum arc [Ref. 2 ], anode spots appear 
under certain operating conditions. These anode spots are 
of practical importance, having been also observed in high 
density plasmas [Ref. 3]. They are typically micron or 
sub-micron sized protrusions on the electrode surface, which 
act to "focus'* the electric field and current stream lines. 
Primaury consideration here is with the anode. 
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B. SBVI2M of PREVIOUS »ORK 
1 . Ph ysic s 

Most work on sheath phenomena is embodied in ”probe” 
theory investigations, where the affects of the electrode 
(probe) locally disturb the quiescent plasma. Such work is 
relevant to this study since the anode is essentially a 
heavily biased probe in contact with the plasma, which is 
quiescent within the sheath region. 

Analyses of electrostatic probes in collision domi- 
nated plasmas were originally done by assuming that a quasi- 
neutral diffusion control led region extended to within one 
mean free path of the probe where it was match ed to the edge 
of the free fall-she ath. No provision was made for a tran- 
sition region between the two regimes [Ref. 4]. 

Shultz and Brown [Ref. 5] studied the transition 
from a collisionless to collision doohinated sheath through 
charged species collection with a spherical probe. Postu- 
lating different physical models for cases ranging from a 
collisionless sheath to a sheath with many collisions, 
semi-em p irical formu las were obtained in each instance. 
Several other authors attempted matching the solution in the 
collision dominated region directly to that of the free-fall 
sheath using variational principles. 
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The previously mentioned techniques can be criticized, 
because in each of these methods something was forced to 
fit. In other words, at least one derivative was disconti- 
nuous at some artificially imposed boundary, thereby raising 
doubt as to the validity of the results. 

The first systematic analysis of probes in collision 
dominated plasmas were carried out by Cohen [Ref. 6] and Su 
and Lam [Hef. 7], who assumed the continuum equations valid 
throughout the plasma, including to the probe surface. As a 
consequence of these analyses, both the sheath and the qua- 
sineutral regions appeared naturally from the diffusion 
equations in the limit as the ra.tio of probe size to Debye 
length approached infinity. The continuum model removed the 
concern about forcing sets of equations to match at the edge 
of an arbitrarily assumed sheath. These authors restricted 
the problem to an isothermal plasma, but as will be shown 
here, the electron temperature may vary significantly in the 
neighborhood of the probe for most discharges of interest, 

Hadbill [Ref, 8] quasilinearized the continuum equa- 
tions, and solved with numerical integration, obtaining 
results that agreed with those of Cohen where comparison 
could be made. 
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In the previous analyses, not only was the isother- 
mal plasma assumed, but the plasma was so slightly ionized 
that only charge- neutral interactions needed to be taken 
into account. However, for highly nonequilibrium plasmas, 
charge-charge collisions play an important role [Refs, 9 and 
10]. The effects of collisions between charged particles 
manifest themselves as volume ionization, volume recombina- 
tion and the dependence of the transport coefficients on 
charged particle density and temperature and, therefore on 
position. Barad [Ref, 11] cites a specific example, 

"If the electrons have a temperature of 2500 ®K, then the 
ionization fraction equal to about 10 ”S will result in 
charge-charge collisions occurring one-tenth as frequently 
as charge-neutral collisions. It is shown... that for 
such a degree of charge-charge interaction, significant 
effects are felt," 

Chung, Talbot and Touryan [fief. 1], in an extensive 
review of probe studies, state that no general solution is 
available for determining charge density and species temper- 
ature for probes small in comparison to the hydrodynamic 
boundary layer thickness. 

Dolson [Ref. 12] investigated the nature and extent 
of the voltage drops in the vicinity of magnetohydrodynamic 
(MHD) non-emitting electrodes, in particular the losses 
attributable to the sheath. tinder his conditions the nonex- 
istence of a one-dimensional sheath solution is shown, and a 
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computer model with two-dimensional, periodic active sites 
representing a flat plate electrode is developed. With this 
model the effects of a magnetic field and Joule heating are 
studied and the results are compared with experimental 
observations reported bj argryoupolos (1973), and Sonju and 
Teno (1974), (Befs. 13, 14]. 

An extension of Dolson's work is reported in 1980 
[2ef. 15], for application to ionizar/sustainer discharges. 
The role of electron pressure is introduced along with net 

V 

ionization and recombination, which proves necessary in 
establishing truely boundary-layer behavior for the charged 
species profiles. Additionally, a modified Newton-Ha phson 
computational procedure, which enables much faster solution 
convergence is formulated. That study demonstrates solu- 
tions with non-active wall densities depressed from the free 
stream density, a result Dolson did not accomplish. Results 
were for 10 volts of electrode potential fall, an increase 
from 5.0 volts reported by Dolson. 

2. Computati onal Bet ho ds 

The fundamental difficulty in the development of 
computational models for the prediction of plasma flow phe- 
nomena lies in the coupling of equations governing both the 
gasdynamic flow and the electromagnetic (EM) fields. Since 

\j 
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the flow fields can affect the EM fields (and vice-versa ) , 
the nonlinear partial differential equations must be solved 
simultaneously and the formulation of the problem becomes 
complicated. In the present study the flow properties are 
decoupled from the the charged species conservation equation 
and Poisson's Equation. 

The approaches used in the solution of plasma flow 
phenomena are quickly reviewed below so as to point out the 
background of the the approach in this study. This problem 
is typical to those dealing with plasma-boundary layer phe- 
nomena. Approaches to solving the general plasma flow prob- 
lems have consisted of the following; 

(a) One dimensional "influence coefficients." 

(b) One-dimensional closed form analytic solutions to 
ordinary differential equations. 

(c) similarity transformation and closed form solutions 
to a modified differential equation in one- or two- 
dimensions . 

(d) One-, two-, or three-dimensional numerical solution 
to the uncoupled steady plasma flow phenomena 
aquations. 

(e) One-, two-, or three-dimensional numerical solution 
to the coupled steady plasma flow phenomena equations 
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using either the magnetic induction equation or the 
elliptic equation governing the electric fields. 

(f) One-, two-, or three-dimensional transient numerical 
solution to the coupled plasma flow phenomena 
equations using either the magnetic induction or the 
elliptic equation. 

Methods (a) through (c) are exact but are restricted 
to a narrow range of problems. Approaches (d) and (e) pro- 
vide rapid convergence to a steady state solution but do not 
consider the transient behavior between the gasdynamic and 
electromagnetic fields, or often the strong coupling between 
the two sets of fields. Method (f) is most desirable since 
it treats a wide range of problems, but suffers from large 
computational times [Ref. 16]. The approach of this study 
is characterized as a method (e) approach. 

C. OBJECTIVES OF THE PRESENT STUDY 

The objective of this study is to demonstrate a numeri- 
cal procedure which can model the physical behavior of a 
conducting fluid, as would be generated in the after-glow of 
an electron beam discharge laser, under the influence of 
electric forces. Calculations are performed in two dimen- 
sions, with a flat plate non-emitting electrode (anode) 
which has periodic active sites. These active sites 
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facilitate current constriction required to satisfy charge 
continuity. Ohm's law, and Poisson's equation. The physical 
model is idealized in Fig. 1. 

The controlling equations are solved for the potential, 
the charged species density distributions, electron tempera- 
ture and current in the vicinity of the electrode. The 
resulting computer programs generate sheath and ambipolar 
regions in a self-consistent fashion, using the same set of 
equations throughout the field. The size of the sheath and 
the voltage drop attributable to its shielding effects are 
investigated. 

This work does not consider those phenomena that lead to 
the arc discharge, which results from thermal instabilities. 
The constriction at the anode is purely that required to 
satisfy continuity. 
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II. PROBLEM FO RMULATIO N 



A. GOVERNING EQUATIONS 

The governing equations that describe the electrode 
region portion of the plasma, consist of Poisson's equation, 
species conservation and electron energy equation, in addi- 
tion to the overall continuity, momentum and energy equa- 
tions of the plasma. These overall plasma equations are 
essentially those of the background gas since the plasma is 
weakly ionized (less than 0.01 %) . 

He have [ Ref . 1 ] ; 

Poisson equation. .. 






( 1 ) 



electron conservation... 






( 2 ) 



ion conservation... 




(3) 
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In gas discharges, the electron energy distribution is 
known to be non-Maxwellian. In order to allow for non-Max- 
wellian distributions and still retain use of a temperature 
concept, an empirical fit to the calculated relation of 
electron temperature as a function of E/n is used in the 
representation of the electron temperature (see Sec. III. A) . 
The concept of a "two-temperature” plasma is preserved but 
no restrictions are made on the energy distribution function 
of the electrons. In addition to the assumption of a func- 
tional dependence of the electron temperature upon E/n only, 
there is no representation of an external magnetic field 
and, in particular the internal or self-generated magnetic 
fields are not considered. The overall gas temperature is 
assumed constant and in equilibrium with the plasma ions. 

The bulk energy equation is assumed to have little bearing 
on this study, since there is negligible gas heating as dis- 
cussed in Section III. A. This effectively uncouples overall 
energy consideration from the problem. 

Tire electrodes serve the purpose of sustaining the S/n 
to pump the laser medium. No account is made of an external 
E-beam or other source of primary ionization. Simply, the 
medium may be considered to exist in an afterglow state. 
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The production term (i.e., the right hand side) in Eqns. 
(2) and (3) has been described utilizing measured and theor- 
etical results of the rate coefficients for Nitrogen ioniza- 
tion and two-, and three-body recombination. The production 
term is expressed as; 

where "i?;" and are respectively the ionization and 

recombination rate coefficients. 

3. BOONDAHI CONDITIONS 

In the customary bounda ry- value problem in fluid mechan- 
ics, a two-dimensional flat plate is assumed. A similar 
two-dimensional Cartesian description is used in the present 
work, which has been found to be the minimum suitable 
description for the problem posed, [Refs. 3,12,15]. No 
one-dimensional solution exists for the frozen flow of 
charges through a collisional sheath, since a Cartesian geo- 
metry in one dimension is devoid of the necessary geometric 
decrease of current density away from the electrode, as 
there would be in a cylindrical or spherical geometry. 

Convection is considered negligible, as can be demons- 
trated by comparing fluid velocity with the drift velocity. 
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Electron and ion densities are determined at the plasma 
edge, knowing the state of the medium. A current balance 
can be used to establish, the species densities resident at 
the active portion of the anode. An electrode voltage drop 
is determined from the properties of the medium; e.g., the 
effective ionization energies of the gas in particular. 
Since current constricts, there must be regions of the flat 
plate that receive no current. The inactive portion of the 
flat plate merits attention, especially in the regions adja- 
cent to the active anode site. Species concentration along 
the inactive portion of the wall can vary, depending on the 
reactivity (catalytic properties) of the wall. There is a 
balance that can be established between the voltage varia- 
tion normal to the wall, and the species concentration at 
the wall together with the species variation normal to the 
wall. A zero concentration of electrons and/or ions at the 
wall is an acceptable boundary-condition; which is to say 
that the wall acts as a perfect sink for the charged spe- 
cies, A detailed kinetic-theory analysis shows that there 
is in fact a lower bound to the electron and ion concentra- 
tions at the wall [Ref. 17]. In any case, there is no cur- 
rent into or out of the inactive wail, (see Fig. 1) . This 
is discussed later in Sect. IV. 
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The remaining portion of the boundary conditions refer 
to the free stream medium, and are easily established. 

C. PARAHETERS OF THE DISCHARGE HEDIDH 

The application of this study is for molecular Nitrogen 
at one amagat, where the gas temperature is 273 <>K. Other 
pertinent parameters follow; 

Slecnric Field: 1.5 to 15 KV/cm. 

Electron Density: to lO**’ m~3. 

Electrode 7oltage Drop: 35 volts. 

The electron energy relation is, as mentioned, a curve fit 
from experimental results [Ref. 18], carried out by R.H. 
Compton and D.J. Sutton (1952). Since cross-section data 
for electrons in Nitrogen discharges are available [Ref. 
19], the electron temperature and diffusion coefficient are 
readily determined, and representable as a function of E/n. 

The ion temperature is assumed to follow that of the 
neutral gas. An electrode potential fall of about 35 volts 
may be assumed since it is the effective ionization energy 
for Nitrogen gas and, therefore, the anode fall [Ref. 20]. 

The cathode is more complicated to describe because of 
the requirement of surface electron emission. The cathode 
fall in cold cathodes is of the order of 250 volts, and is 
considerably less in thermionic cathodes. A proper 
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III. PHO BLE H SOLUTION 



A. PSOBLEH SIMPLIFICATION (PHYSICS) 

There are a number of simplifications that are incorpo- 
rated in the problem description. Namely; 

• Steady state conditions (within the pulse time). 

• No magnetic fields. 

• Negligible convection. 

• Negligible Joule heating, T^> = T^= constant. 

However, specific account is made of; 

• lonization/recombina tion. 

• Electron pressure. 

These are discussed in part in this section. 

1 . Conv ecti on 

The aabipolar region is the transition region from 
the sheath to the undisturbed plasma and can perhaps span 
the boundary layer. This fact requires some consideration 
of the presence of convective effects. 

From the ion species conservation equation (Eg. 3) , 
for example, we know that in the ambipolar region the last 
term of the left hand side is small compared to electric 
conduction. Now we need only to compare the convection with 
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the conduction, or simply the fluid velocity with the drift 
velocity. 

-yi 

V = ^ 0.1 to 1.0 m/sec 

where the Reynolds number (R^) is 10^ to 10‘ and the free 
stream velocity is 10 to 100 m/s. Now the drift velocity 
may be approximated by 

Vj as •= 1 0 to 100 m/s 

for E = 105 to 10* V/m and the mobility JU. = 10~* 
a2/s/volt. Clearly, in the presence of a relatively strong 
intarelectrode field, the contribution of ion convection due 
to a cross flow should be negligible. The geometrical 
orientation of the discharge with respect to the primary 
flow can thus become an important consideration. 

2. D ensit y B oundary La ver Forma tion Time 

The charged-species density profiles change in a 
more complicated manner than either the voltage potential or 
the electric field. The magnitude of the ion and electron 
densities may change appreciably within the sheath; more- 
over, fractional analysis is risky because inflections are 
present in the profiles, k discussion pertaining to the 
stabiliry of the ambipolar region and its boundary layer 
nature is presented in detail by Biblarz, et al. , [Ref. 15]. 
However, Appendix A shows that the characteristic length of 
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the ambipolar region is of the order of the characteristic 
sheath length. Letting L=5x10~® • meters and the drift veloc- 
ity as before, v^ = 10 to 100 m/s, a characteristic time for 
the formation of the density boudary layer would be 0.5 to 
5.0 microseconds. It is reasonable that the ambipolar dif- 
fusion to the wall can establish itself well within the 50 
microsecond pulse period of the electrodes. 

3* Si g nificance Of Joule H eati ng 

For the present problem, joule heating can be shown 
not to be significant. In the vicinity of the active site, 
the average energy deposition rate per unit volume is shown 
to be approximately 4x10® watts/m^. The temperature of the 
Nitrogen gas at one amagat would be raised approximately 1°K 
in 2.6x10-* sec. Allowing a 50x10-* sec pulse period would 
raise the temperature of that region by approximately 20®K, 
or 7%. 

Assume that the heat generated at the anode spot is 
carried away by thermal diffussion (conduction). Then for 
Nitrogen gas ()c=0.0267 H/m-®K, f =1.138 Icg/m^, and Cp=1043 
J/kg-®K) the thermal diffusivity is; 

0 (s -^ = 2.25x10-3 m2/sec 

A typical diffussion time would be (for, Aj =4.4x10-® m) ; 

0.85 millisec 
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The "flush time" is the time required for clearing 
the lasing cavity of residue gases in preparation for the 
next pulse, and can be estimated as the gas flow rate 
divided by the flow length. Assuming v = 100 m/sec and L = 
10 cm, the flush time would be; 

FT = 0. 1/100 = 1 fflillisec 

Other factors considered, the inter-pulse period (time bet- 
ween pulses) would be larger, but even in a one millisecond 
period, there is just enough time to cool the anode spot. 

The analysis so far has assumed a flow velocity at 
the surface to be the same as the bulk flow in the laser 
cavity. However a boundary layer exists — probably turbu- 
lent, with a laminar sub-layer adjacent to the electrode 
for ; 

Rg = X^/^ = (0.1) (100)/( 15.53x10-6) = 6.44x105 

The velocity boundary layer thickness would be 2.55 millime- 
ters. Assuming a 1/7-th power profile, the velocity at a 
sheath length from the electrode would be: 

V = V(\/5/^= 100 (4.4x 10-5/2.55x10-3)1/^ = 55 m/sec 
A characteristic local flush time would be approximately 2 
milliseconds, and much higher for a laminar velocity pro- 
file. These flush times can be quite large as the require- 
ment to cool the Joule heated gas near the anode spot is 
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increased. The design flush time based on other factors may 
infact be insufficient to clear enough of the Joule heated 
gas from the vicinity of the anode spot, enabling a build up 
of heated gas, which could lead to local enhanced ionization 
of the gas. 

4. E lectron Temp e rat ure 

The electron temperature (T^ ) is assumed to be a 
linear function of E/n as presented in Pig. 2, for the 
region of interest, namely E/n = 0.5 to 5 . 0 x 10 “ 2 o v-m^. For 
the range of S/n considered, it is adequate to take the 
electron temperature as; 

Te = [A«log(E/n) + B]*T<, (5) 

where; 

A = 12. 1 

B = 38. 5 

E/n is normalized to 10 * 2 o 
The diffusion coefficients were determined from 
Brown [Ref. 18]. The drift velocities were obtained from 
plotted results, and since /^=vj /E, and Pj=/ 4 ^kTj/e, (Ein- 
stein's relation), the diffusion coefficients are; 

Dg = 0. 0665 m2/sec 
Dj = 8. 0x10~* m2/sec 
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B. 



HORKIUG SQOATIONS 



Since convection is negligible, and assuming that the 
sustainer operates in the afterglow of an electron beam, the 
governing equations, namely Eqs. 1,2 and 3, can be written 
as follows; 




( 6 ) 



~y*[Vn^+ /Q - n^V<?^/d] =C2.ng 



(7) 



-7*[7n- + 






c,n^ 



( 8 ) 



where; 



A. 




e^/kT^ 


e = 







c, = eiji/kTj^ Cj^= Aj^n«^/D^ Cj= Cj^D^/D^ 

The equations may now be applied to flat plate formula- 
tion as depicted in Fig. 3. It is convenient to drop the 
•w\rt for simplicity, but it is to be understood that all 
variables have been suitably no n-dimensionalized. 

We have Poisson's equation; 







c^jrIc■'n 






(9) 
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electron conservation equation; 
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( 10 ) 



ion conservation; 









(11) 



C. NOMERICAL METHODOLOGY 

The nonlinear terms creat problems in the computer anal- 
ysis of the electron and ion species conservation equations, 
(Bqs. 10,11). The Jacobi method includes all nonlinear 
terms on the "right hand side"; i.e., external to the coef- 
ficient matrix. Convergence to a solution is possible if 
these non-linear terms change slowly over each iteration. 

It was found by Dolson [Ref. 12] that the Jacobi method was 
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in fact unstable for the present set of equations and 
conditions. As a consequence of the failure of the Jacobi 
method, a quasi-Jacobi method wcis implemented, in which an 
estimate for each of the solution variables n^^) was 

computed. When the product of two variables is encountered, 
one variable is treated as a constant coefficient for each 
iteration. This means that the non-linear terms are 
retained in the coefficient matrix. The "constant” coeffi- 
cients are updated after every iteration, thus altering the 
coefficient matrix. The conventional Jacobi method was 
found to converge only for low values of voltage potential 
at the electrode, whereas the quasi-Jacobi procedure pro- 
vided converged solutions for approximately 5 volts. A dis- 
cussion of these methods is provided in Appendix C as back- 
ground for the development of the modified Jlewton-Raphson 
(MMB) method, which is used in the solution to Sqs. 1 
through 4. 

The nonlinear coupled partial differential equations 
(PDB) , Eqs. 9,10 and 11 are linearized to provide a set of 
linear algebraic equations, which must hold throughout the 
domain described as follows. The domain as illustrated in 
Pig. 3, is divided up into equal sided sub-areas whose ver- 
tices are the computational nodes at which the linearized 
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algebraic equations must hold. The wall is represented by 
seven (7) nodes, with the leftmost node being the location 
of the active site. The "y" direction extends from this 
active site along the wall. The "x” dimension perpendicular 
to the wall is represented by 25 nodes starting from each of 
the wall nodes. The solution to the PDE's is represented by 
the vector z(x,y) , or on the computational mesh as z (i, j) . 

The boundary conditions have been discussed previously 
and are presented in Pig. 3. The boundary condition that 
there is no net current into the nonactive wall (^=0) is 
satisfied by the electron contribution to current only. 

This is justified since it can be shown that the electron 
contribution is approximately three orders of magnitude 
greater than the ion contribution to the total current. 
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Boundary Conditions for Computational Method. 
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IV. EBSOLrS MD CONCLUSIONS 
A. GENERAL DISCUSSION 

Computer solutions are reported here for three cases 
presented in Table I. The active site voltage was 35 volts 
in each of the three cases. Included in this section are 

TABLE I: Computer Solution Conditions 

CASE ng(wall) n; (wall) n^ (anode) n^ (anode) 

I 1.0 0.0 J-matched 0.0 

II "float” 0.0 not J-matched 0.0 

III "float" "float" not J-matched 0.0 



also some pertinent results obtained at 10 volts to illus- 
trate the effects the net production term has on the conver 
gence to a boundary layer behavior for the species density 
and free stream electric field. 

The solutions are depicted as obligue presentations of 
the potential, species density, space charge, net produc- 
tion, electric field and Joule heating distributions over 
the "x-y" coordinates defined in the problem formulation. 
Also included are two-dimensional plots of various terms 
along a "cut" or line perpendicular to the electrode wall. 
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Suart.ing at the active site, with some starting at the wall 
along the inter-nodal line of symmetry (see Fig. 3) . The 
solutions are organized as Case I,ri and III results, will 
be presented as Pigs. 12,13 and 14 respectively at the end 
of this section. The discussion of these results proceeds 
with additional figures which are interspersed in the text 
appropriately. 

Table I presents various electron and ion boundary con- 
ditions along the wall some of which are "floated." By this 
term is meant that upon the completion of each numerical 

ti 

iteration, the electron (and/or ion) density at the nonac- 
tive portion of the wall is re-established by assigning it 
the value of a near neighbor away from the wall. This pro- 
cess is illustrated in Fig. 4. 






Figure 4 

Description of "floating the Wall Boundary Condition. 
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The usual elliptic boundary value problem has well 
possed boundary conditions from a priori information. How- 
ever, the nature of the problem in this study includes very 
little information as to the correct species wall density. 

It seems reasonable to set these densities to zero along the 
inactive portion of the electrode as previously mentioned in 
Sec. II. B. In this study the wall boundary conditions are 
"lossened" to facilitate "naturally occuring" and unpresum- 
ing species density profiles. The nature of the problem in 
effect, produces a self-consistent boundary condition. 

B. EFFECTS OF NET PRODUCTION 

Earlier work [Ref. 15] reported limited results for a 10 
volt active node site, solutions for cases with and without 
the inclusion of the net production term have profound 
influences on the appearance of the charged species density 
profiles. Figures 5,6 and 7 illustrate the before and after 
effects of including the net production term. In both cases 
the wall boundary values for the species density was ini- 
tialized at the free-stream value (i.e., 1.0), then allowed 

to "float" as described previously. The addition of the net 
production term brings about the boundary layer behavior of 
the densities of the electrons and ions is illustrated in 
Fig. 7. Efforts to attain converged solutions with the net 
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pcoductz.on tern set to zero proved fruitless for active site 
potentials exceeding 10 volts. 




Figure 5 

(a) Slectron and (b) Ion density perspectives with 
the net production tarn set to zero. (10 volts) 
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Figure 6 

(a) 3lectron and (b) ion density perspectives with 
the net production tern coupled. (10 volts) 
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Figure 7 

Influence of the net production term on the 
density profiles along a line perpendicular to the wall, 
extending from the active node site, (10 volts). 
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C. DISCOSSION OP RESULTS 



1 . Voltag e Potenti al 

The voltage potential of the active node site ("the 
anode") in each of the cases presented was* 35 volts. Fig- 
ures 12-a, 13-a and 14-a will clearly illustrate sharp vol- 
tage spikes at the anode, and the "Laplacian" nature of the 
field far from the anode. The potential along the inactive 
wall varied from 12 volts at the farthest position from the 
anode to 17 volts just next to the anode. Similar profiles 
of the electric potential along the wall for the other cases 
are presented in Fig. 8. Figures 12-i, 13-i and 14-i will 
show the voltage potential variation along a line normal to 
the wall extending from the anode site. 

Figure 9 illustrates a two-dimensional cut of the 
voltage potential along the inter-nodal symmetry line. Here 
it can be seen that the variation from the wall is shallow. 
It is important to point out that the electron temperature 
is directly related to the gradient of the voltage potential 
through Eg. 5. This fact has certain consequences which are 
discussed in the next section. 
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Figure 8 

Voltage potential along the wall boundary (electrode) 
for CASES I, II, and III. 
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Voltage potential variation normal to the wall, 
along the inter-nodal symmetry line, (CASE I) 
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2* e lectro n EnergY 

The electron energy will be presented in series -b 
of Figs. 12,13 and 14, as oblique perspectives. The anode 
site temperature for Cases I, li and III are 2.1, 2.1 and 
1.9 e7 respectively. The description of the electron temp- 
erature is that of a linear function of the logarithm of the 
electric field. Any portion of the potential field that is 
’’level” (grad <P =0) will necessarily generate a low or even 
negative (see Fig. 14-b) electron temperature. According to 
Fig. 2 and the discussion in Sec II. A. 4, a negative tempera- 
ture would occur below an B/n of 0.04x10-2o v-cm^. This 
problem can be removed by setting a minimum value to the 
non-dimensional temperature field, nominally 0.13, corres- 
ponding to the minimum at the low end of E/n in Fig. 2. 

The parameters of the discharge (Sec. I.C) require 
an electric field of 2.7x10s v/m (or, log£E]=5.43) at the 
free stream boundary, which corresponds to an electron temp- 
erature of 0.905 eV (^=1.) . The temperature fields in all 
cases presented do not quite meet the free stream boundary 
requirements. The electric fields are insufficient to ren- 
der proper results to the electron temperature profiles. 

Thus the free stream temperature boundary requirements are 
"forced” by slightly inflating the influence of the electric 
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fisii. The effect of forcing the temperature field in this 
way can also be looked at as "sliding*' the linear 
representation of the temperature curve to higher values at 
constant S/n, (see Fig. 2). Figure 10 illustrates a case 
where a low electric field, here 0.8 kv/cm, results in a low 
electron temperature projecting toward the free stream 
(d=0.62). The low temperature is adjusted to meet the 
desired free stream value (5=1.0) by multiplying the elec- 
tric field by a factor of 3.40. Figure 10 also illustrates 
an alternate method of meeting the free stream conditions by 
truncating the variation of temperature with a minimum 
value, namely the free stream boundary condition. 

Both of the "forcing" techniques are used, but for 
different reasons. As long as the adjustment factor does 
not exceed an order of magnitude change in the electric 
field, then it is not objectionable in light of the simpli- 
fied representation of the electron temperature. The 
adjustment requirements of the three cases reported in this 
study are between 3.0 and 5.0, with Case I using a factor of 
10 to see the effects on the Joule heating. In Case I the 
non-dimensional temperature is artificially high at 1.4 eV 
at the free stream. 
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Figure 10 



Electron temperature variation normal to the wall, along the 
nodal symmetry line; illustrating two methods of , ad justment 
to "force" the free stream temperature conditions. 
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The electric field near the inter-nodal symnetry points 
of the wall are characteristically 2 orders of magnitude 
less than that of the free stream conditions. The resulting 
very low, even negative temperatures in Case III, Pig. 14-b, 
is *»patched" using the truncation method previosly dis- 
cussed. Whether these low temperatures and the associated 
steep derivatives are physically reasonable, merits further 
investigation. 

The temperature fields generated and presented as 
Figs. 12, 13, 14 series -b are the final result depictions. 
However, the temp^ature fields used in the solution process 
were adjusted and truncated as previously discussed to main- 
tain solution convergence. An example of a "corrected” 
temperature field is presented in Fig. 14-c. 

The electric field and Joule heating distributions 
are generated from the primary voltage and electron density 
solutions, much in the manner as the temperature distribu- 
tion was determined. In fact, the Electric field and the 
Joule heating distributions are similar in appearance to the 
temperature distributions. These distributions play no 
direct part in the solution scheme and thus being ancillary 
in nature were not "forced” to meet any a priori conditions. 
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The Joule heating tern can be written as; 

= Ci>en^^icr.^/A^2].C^ -7ng]*V^ 

and thus is a rather straight-forward function of <fi, n^, 
and a,-. The ion contribution is shown to be negligible, 
since . The Joule heat at the anode site is of order 

10® and is 100 tines that of the free stream values. Joule 
heating effects are negligible in the overall problem, as 
discussed in Sec. II. A. 3, obviating the reguirenent for an 
energy equation. 

4, Specie s Density Dis tri bution 

The solutions for the charged species density dis- 
tributions are presented as oblique perspectives in Figs. 
12,13 and 14, series -e (electron density); series -f (ion 
density) ; and series -g (space, charge density) . Addition- 
ally the variation of these distributions normal to the wall 
along the active node symmetry line are presented in the 
series -j and -k figures. The oblique perspectives of the 
densities are reversed for clarity. 

The electron and ion densities at the nonactive wall 
are a subject of interest in the literature, to which little 
information has been gathered on a firm basis. The results 
of this study suggest a slightly positive space charge at 
the inactive wall. Scrutiny of Case III results reveals, a 
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2.2x10‘* du“ 3 positive space charge number density, (or 3.5 
milli-coulombs) . Recall that Case III involved ''floating” 
both of the charged species densities at the wall. In fact, 
close inspection of Case I and II reveals a positive space 
of similar magnitude but displaced slightly from the wall. 
The wall in these cases was set with a negative space charge 
as the boundary condition. The positive space charge along 
the wall is most predominant in the vicinity of the inter- 
nodal symmetry line, and whether the low temperature values 
in this region have significance, remains to be determined. 

The active node site (anode) electron density in 
Case I (2.4x10^^ m-3) was determined so as to match the free 
stream current to the current constricting at the anode. In 
the results of Cases II and III the anode electron densities 
are chosen so as to provide a smooth variation extending 
outward from the active site. 

5. l oniza t ion/Recomb in atio n 

The net ionization term is utilized in ail cases; 
and was found necessary for (1) establishment of the charged 
particle density distributions, and (2) solution convergence 
for sheath lengths spanning more than one computational node 
spacing. The net production term is depicted in Figs. 12,13 
and 14 series -h. There is net ionization throughout the 
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solution domain. However as Case III demonstrates, there is 
a factor of 5 orders of magnitude greater ionization within 
the sheath region surrounding the anode site than at the 
free-stream edge. Cases I and II show ionization along the 
nonactive wall, simply because the ion density is zero, eli- 
minating significant recombination, (see Eg. 4) . 

6 • Current D escr iption 

Figure 11 is typical of the constricting nature of 
the current stream lines. The potential, electron density 
and the electron temperature are used to determine the ”x-y” 
components of the current density. Using a “nearest neigh- 
bor" influence principal, a local derivative can be deter- 
mined; from which a streamline can be "marched-out" in small 
increments. The profiles of all three cases are similar, 
the fact being that the current matching discussed previ- 
ously will not affect the current streamline profiles. Fig- 
ure 11 clearly demonstrates the satisfaction of the current 
constriction reguirement necessary to permit solutions to 
the system of eguations. 

The current into the anode site in Case I is matched 
to the current entering at the free-stream boundary. This 
specifies the electron density at the anode site. However, 
in Cases II and III the current at the anode was only three 
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tines the free-stream currant. The current matching problem 
does not deserve serious complaint, since the problem 
description is simplified tc two-dimensions and the current 
constriction is only reguired to satisfy continuity. 




Figure 11 
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Figure 12-a 



Figure 12-b' 
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Figure 12-c 



Figure 12-d 



Oblique perspective of 
log Electric Field (CASE I) 



Oblique perspective of 
log Joule Heat (CASE I) 
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Figure 12-e 

Oblique perspective of 
alectron uensity (CASE I) 
(reversed perspective) 



Figure 12- f 

Oblique perspective of 
Ion Density (CASE I) 
(reversed perspective) 





Figure 12-h 

Dblique 'perspective of 
Net Production (CASE I) 



Figure 12-g 

Oblique perspective of 
Space Charge (CASE I) 
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Figure 12-i 

Presentation of the Voltage Potential 
(CASE I) along a line normal to the. wall, 
extending rrom the active node site. 
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Figure 12-j 

Presentation of Species Density 
(CASE I) along a line normal to the wail, 
extending from the active node site. 



57 



SPRCE CHRRGE DENSITY 



1.0 



0.8 



0.6 




DISTRNCE (MM) 



» 



0. 15 



Figure 12-k 

PresentatioQ of the negative Space Charge Densit 
(CASE I) along a line normal to the, wall, 
extending from the active node site. 
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Figure 12-1 

Presentation of the Electric Field 
(CASE I) along a line normal to the wall, 
extending rrom the active node site, 
(logarithm scaled) 
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Figure 12-m 

Presentation of the Joule Heating 
(CASE I) along a line normal to the. wall, 
extending from the active node site, 
(logarithm scaled) 
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Pigura 12-n 

presentation of the Species Densities 
(CASS I) on a line normal to the wall., 
extending along the inter-nodal symmetry line. 
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Figure 13-a Figure 13-b‘ 

Oblique perspective of Oblique perspective of 

Potential Field (CASE II) Temperature field (CASE II) 




Figure 1 3-c 



Figure 13-d 



Oolique perspective of 
log Electric Field (CASE II) 



Oblique perspective of 
log Joule Heat (CASE II) 
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Figure 1 3- e 



Figure 13-f 



Oblique perspective of 
Electron Density (CA5S II) 
(reversed perspective) 



Oblique perspective of 
Ion Density (CASE II) 
(reversed perspective) 




Figure 1 3-g 

Oblique oerspective of 
Space Charge (CASE II) 



63 



Figure 13-h 

Oblique perspective of 
Net Production (CASE II) 
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Figure 13-i 

Presentation of the Voltage Potential 
(CASE II) along a Ixne normal to the wall, 
extending from the active node site. 
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Figure 13-j 

presentation of Species Density 
(CASE II) along a line normal to the wall, 
extending from the active node site. 
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Figure 13-k 

Presentation of the negative Space Charge Density 
(CASE II) along a line normal to the wall, 
extending from the active node site. 
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Figure 13-1 

Presentation of the Electric ^-i 

(CASF. HI alona a line normal to th^ wall» 
^extending from the active node site, 
(logarithm scaled) 
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Figure 13-in 

Presentation of the Joule Heating 
(CASE II) along a line normal to the wall, 
extending from the active node site, 
(logarithm scaled) 
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Figure 13-n 

presentation of the Species Densities 
(CASE II) on a Ixne normal to the wail,, 
extending along the inter-nodal symmetry line. 



69 




Figure 14-a 



Figure 14-b 



Obligue perspective of 
Potential Piela (CASE III) 



Dbligue perspective of 
Temperature field (CASE III) 




Figure 14-c Figure 14-d 

Oblique perspective of adjusted obligue perspective of 
Temperature Field (CASE III) Joule Heat {CASE III) 
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Figure 14-e 

Oblique perspective of 
Electron Density (CASE III) 
(reversed perspective) 



Figure 14-f 

Oblique perspective of 
Ion Density ^CASE III) 
(reversed perspective) 






Figure 14-g 

Oblique perspective of 
Space Charge (CASE III) 




Figure 14-h 

Dblique perspective of 
Met Production (CASE III) 
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Figure lU-i 



Preseatation of 
(CASE Illi along a 
extending from 



the Voltage Potential 
line norm to the. wall, 
the active node site. 
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Figure 14-j 

Presentation of Species Density 
(CASE nil along a line normal to the wall, 
extending from the active node site. 
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Figure 14-k 



Presentation of the negative Space Charge Density 
(CASE III) along a line normal to the wall, 
extending from the active node site. 
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Figure 1U-1 

presentation of the Electric Field 
(CASE III) along a line normal to the wall, 
exten^ng from the active node site. 

(logarithm scaled) 
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Figure 14-« 

Presentation of the Joule Heating 
(CASE III) along a line normal to the wall, 
extending from the active node site, 
(logarithm scaled) 
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Figure 14-n 

presentation of the Species Densities 
(CASE III) on a line normal to the wall, 
extending along the inter-nodal symmetry line 
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V. C ONCLUS IONS AND RECOMM EN D ATION S 



A. REVIEW OF ASSUMPTIONS 

In the development of this work several assumptions were 
made which the solution results later verified. The assump- 
tion of negligible convection effects is straight forward. 
The charged particle density boundary layer evolution time 
turns out to be small in relation to the pulse period; thus 
enabling consideration of the problem as a steady-state aft- 
er-glow problem. The fact that the Joule heating of the gas 
is shown negligible, frees the formulation from an elaborate 
energy description. since the electron temperature is a 
simple function of E/n, the problem can be further simpli- 
fied. The solution proceeds with nominal values of the par- 
ameters of the discharge as discussed in Sec. II. B. 

Finally, the importance placed on the net production term is 
verified by the profound changes in the results including 
the fact that ”n^" is necessary for the establishment of the 
charged-species boundary-la yer behaviour. The analysis of 
App. 3 shows that an assumption of "frozen flow" would in 
fact be premature. 
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B. PHYSICAL CONCLUSIONS 



The present model of the sheath and ambipolar region 
evolves from the assumptions of steady state, uniform ion 
temperature and constant diffusion coefficients. A summary 
of important results is presented in Table II. 



TABLE II; SUMMARY OF IMPORTANT RESULTS 



SUMMARIZED 

RESULTS 


CASE 

I 


CASE 

II 


CASE 

III 


(units) 


sheath length 


3.3 


3.5 


3.7 


10-5 m 


ambipolar length 


4.5 


3.5 


3.7 


10-5 31 


E^(f ree- stream) 


7.9 


3.9 


6.3 


10 ♦ V/m 


(anode) 


1.9 


2.5 


2.2 


10* V/ra 



The results show that a sheath develops as anticipated 
from theory and the proposed assumptions. The sheath char- 
acteristic length is consistent with that predicted by Deby 
shielding. Debye shielding of a flat plate predicts that 
the sheath length is 4.4x10“5 meters. For a poinr (or 

spherical) shielding the sheath length is 3.1x10“5 meters. 
The present study involves the shielding of a point voltage 
source situated on a surface charged plate. Since both two 
and three-dimensional effects are inherent in the formula- 
tion, the resulting sheath lengths are quite reasonable, 

u 
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As Table II shows, the size of the ambipolar region is of 
the same order as the sheath size. Here the ambipolar 
length is defined as the length extending from the edge of 
the sheath to the point at which 99 % of the free-stream den- 
sity is recovered. Finally, the sheath can be self-generat- 
ing from a consistent set of equations without an overall 
energy equation. 

The resulting electric fields at the free-stream boun- 
dary do not quite attain the desired 2.7x10® V/m presented 
in the problem formulation. These results fall short by the 
same factors required for the adjustment of the electron 
temperature field to meet the stated free-stream boundary 
requirements. It must be pointed out that the electric 
field is the required boundary condition. However, this 
boundary condition could be made to yield no reasonable 
results, and so the chosen boundary condition is = 0. 

The resulting electric fields at the active site in each 
of Cases I, II, or III was of two orders of magnitude higher 
than that of the free-stream. These electric fields are 
below that needed for breakdown of Nitrogen at one atmo- 
sphere, which is determined as, [Ref. 18]; 

3^ = (E/p) p = (350 V/cm-torr) x(760 torr) = 2.7x10' V/m 
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The ragairement of strictly matched curtents produces 
steep derivatives and inflections of the density profiles in 
the vicinity of the anode. A noticeable "spike" is evident 
in Fig. 12- j for example. Case II and III anode electron 
densities were established on more aesthetic grounds, where 
smoother appearing density profiles were generated. In 
these cases the anode current exceeded the free stream cur- 
rent by only a factor of 3. 

The nonactive regions of the electrode (wall) far from 
the active site reveal a small positive space charge den- 
sity amounting to approximately 2/4 of the free-stream den- 
sity. The negative space charged sheath is 46% to 56% of 
the free-stream density. 

Previous work [Ref. 12] of a similar problem reported no 
reasonablo results for the case of a catalytic wall 
(n^=n^*=0). The results presented in this study demonstrate 
wall conditions of n^- = 0. and n^ < 1 . Here, a noticeable 
positive space charge coexists with the large negative space 
charged sheath, without large inflections in either the vol- 
tage potential or densities. 

The Case III study has both electron and ion wall condi- 
tions floated. The shape of these density distributions is 
guite symmetrical about the anode site in all directions. 
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One needs little imagination to visualize that the density 
profiles can be "wrapped" around the anode completely, in 
effect treating the anode as a whisker or protrusion into 
the plasma. However, the imagination must be tempered by 
the fact that it can only be a two-dimensional protrusion 
(e.g. , a "blade") . 

The effect of the net production term has been previ- 
ously discossed as one of the most important results of this 
study. The term is necessary for inclusion in the problem 
in order (1) to attain wall electron density values of less 
than the free-stream, (2) to produce the boundary-layer 
behavior of the species densities, and (3) to obtain con- 
verged solutions with voltages greater than 10 volts or com- 
putational node spacings of any order less than that of the 
sheath length. 

C. NUMERICAL COHCLOSIONS 

The Modified Newton-Raphson (HNR) method, presented to 
solve the particular set of coupled non-linear partial dif- 
ferential equations in this study, is a general method that 
can be applied to more sophisticated systems. The defi- 
ciency of the guasi-Jacobian method (discussed in Sec. III.C 
and App. C.) in dealing with non-linear terms is circum- 
vented by an alternative linearizing scheme. Convergence is 
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rapid; however, the starting solution must not exceed 10 
volts, must have zero net production, and must have a cons- 
tant electron temperature during the first 2-3 iterations. 
After the solutions take their initial form at ten volts, 
the potential field is "amplified” and reinstated as a new 
starting solution. The final stages of the converging solu- 
tion are carried out under the full influences of (1) float- 
ing wall densities, (2) variable electron temperature, (3) 
net production and (4) an active node site electron density 
being solved for through current matching. 

The profiles of the space charge and net production dis- 
tributions demonstrate erratic behavior at the first few 
computational nodes from the active site. This is believed 
to be caused by the linearization of the governing non-li- 
near equations and has the appearance of truncation error 
when observed on an iteration to iteration basis. However, 
the remaining regions of the profiles are smooth. 

D. HSCOHMSNDATIONS 

This work succeeded in solving the problem in a two-di- 
mensional environment. This fact necessarily limits the 
quantitative results since the current density decreases in 
only two dimensions. With the computational power of recent 
computers, it should prove worthwhile to make the next 
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■r . 









advance as a three-dimensional problem rather than elaborate 
in the two-dimensional case. A three-dimensional scheme 
would render more realistic results of a quantitative nature 
because active spots are more likely to exist than active 
lines and because 3-D geometries require less voltage for a 
given current. 



Specific recommendations other than the above include; 



• Account for the ion contribution to the current 
boundary condition into the wall. 

• Investigate the physical implications of the low 
electron temperature near the wall far from the 
anode site. 

• Investigate the physical implications of the elec- 
tron density "spikes” prevalent near the anode sits. 

• Make the wall temperature a valid boundary condi- 
tion. 

• Use a free-stream boundary condition of 

=2.7x105 V/m, rather than zero potential. This 
would introduce desired results for pumping lasers. 

• Include a more sophisticated electron temperature 
model, or verify the present one. 

• Increase the computational mesh density. 

• Study the general application limits of the MNR 
method. 

• Investigate "n " in the role of making one-dimen- 
sional solutions possible (i. e. , alter the reacting 
rates) . 

• Investigate whether the low electron temperatures 
and associated steep derivatives located along the 
wall between active sites have any physical implica- 
tions, or result merely from inadequate problem for- 
mulation. 

• Consider the electron-beam ionization in the solu- 
tion to produce more realistic effects. 
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APPENDIX A 



CHARACTERISTIC PARAMETERS 

A. SHEATH LENGTH 

The extent of the sheath is one of the most important 
charcteristic lengths in this study. It is within the 
sheath that most of the potential drop occurs for a short 
discharge. Fortunately, the sheath length can be estimated 
rather easily from Poisson’s equation, (Eg. A1) 








Let; 

. ^(41 



^5 - ^s/nec» 



(A1) 



(A2) 



where is the characteristic sheath length; the sub- 
script '*o” indicates the xalue at the electrode and ” " the 

value at the undisturbed plasma. After gathering the dimen- 
sional quantities on one side. Eg. A2 becomes; 




The characteristic sheath length defined in this 

fashion, is used to non-dimensionalize the “del” operator in 
all directions. That is, independent of the dimensionality 
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of the problem, the sheath is effective over the character- 
istic length just defined. Now an estimate of the size of 
the sheath can made by setting the grouping of dimensional 
quantities in Eg. K3 to be of the order one, namely; 




where ” is the Debye length in the undisturbed plasma. 

Anticipated values are now established for the sheath 
length. Take for example a 35 volt anode potential, and a 
charge species density of 10^^ to ra“3. These parame- 

ters define a sheath length of 1.4x10“^ to 1.4x10-* meters. 

As can be seen, a ten-fold increase in voltage would 
change " Xj” a factor of 3, rendering the nondimensionali- 
zation of the problem relatively insensitive to the chosen 
value of the voltage potential. That is to say the value of 
the sheath characteristic length as defined may be consid- 
ered a resonable estimate for the charge species densities 
examplif iad. 

The important implication here is that extreme cafe must 
be devoted to electrode surface preparation to preclude 
whiskers or flaws of a size comparable to These pro- 

trusions would produce a focusing of the electric potential 
to create extremely high electric fields, which become sites 
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for initiation of breakdown. Protrusion growth has been 
observed to fora in a time on the order of milliseconds. 
Many conditions are favorable to protrusion growth, depen- 
dent on the pre-breakdown and breakdown current profile, 
electrode composition (especially in the case of a ther- 
mionic cathode) , surface purity and surface preparation 
[Ref. 21]. Some possibilities arra; 



B. DIFFUSION LENGTH 

If diffusion is a significant factor in this study then 
its contribution to the current would be; 



or, since > Dj , and using Einstein's relation, or 
D/^«=kT^/e; 



making the equation unit free (through fractional analysis) 



rieid due to micro-particie nomoardment. 

• sticking of microparticles. 

• Electron beam induced protrusion growth. 






e l>t <5 x' 



JL _ <^f\f 






(A5) 
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If the right side of Eg. A5 is of the order one, (the nondi- 
mensional density gradient is a maximum of 1) then ”L" may 
be estimated as; 

L eD^n^/J « 3.7x10”* m 

for the values; 

ne = 10‘« m-3 

= 0.0665 m2/(y«sec) 

= J4fe»n^*xi^»{E/n)ss 2.88x103 A-m“3 
Since L, then any role attributed to diffusion 

should be visible within the sheath. The E/n within the 

thin sheath is of course greater than the E/n at the free 

stream (IO 20 v«m2), which decreases the local value for "L" 

still further. 

Another approach is to assume a one-dimensional solution 
to the charged species conservation equations, which can 
only be done outside of the sheath, i.e., in the ambipolar 
region. The ion conservation equation may be written; 

n = 7.J = - 07n ] (A 6) 

Since it can be shown (next section) that the production 
term can be approximated as; 

n jtP^*n -bCn 2 

then for the stationary condition of no net production; 

2^/ =<X.'S 

Sondimensionalizing Sq. A6, it can be shown that the 
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characteristic diffusion length is; 



L = \fO/pfN = 2.x10“* meters 

which is of the same order of magnitude as the characteris— 
tic sheath length. 

C. lONIZATION/RECOMBIHAnON LENGTH 

Analysis of the production term (Eg. 4) reveals that for 
Nitrogen the three-body (electron) collisional recombination 
rate is small in relation to to the other terms, (see 

App B.) . In any case the net production term is; 

In the steady state the net production term is zero, and 
as is done in Appendix B, the relation becomes; 

n^ = (<*',+ (y^n*) *(n^- n^) *n^ 

M 

substituting the value of 0.17x10-ii m^/sec, and (ir - ) cc 

1Q1®, we can estimate the characteristic production time to 
be Tp = 5.88x10-' sec. 

To establish a meaningful characteristic length, the 
average thermal velocity will be used since it represents 
the principal transport velocity of the medium. 

= = 5.88x1 0-'»3. 2x102 = 1,88x10-* m 

The characteristic length for net production is a factor 
of four larger than the sheath length ( 4.4x10~® m) . The 
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coaclasioa which can be made from this (fractional) analysis 
is that the net production is significant within the overall 
region of interest, including the sheath itself. 

D. FROZEN PLOW 

Hhen a process relaxation time is long compared 

with the time that a fluid particle spends in a given 
region, the flow is said to be "frozen" in this region. By 
this term one means that a reaction process may be 
neglected, as far as influencing the flow is concerned, and 
the species involved treated as chemically inert, [Ref. 22]. 
In the opposite limit, where "T)." is very small relative to 
the characteristic flow time, the flow is said to be in 
"equilibrium." For the recombination rates given in this 
section, "7J." ~ 6.xl0”^ and "^" = 4.0x10~^ seconds. Here 
flow time is defined as the sheath characteristic length 
divided by the drift velocty. Thus the flow is not "frozen" 
as far as chemical reaction, but we must consider that the 
electron and ion species are not in thermal equlibrium. 
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APPENDIX 3 



DEVELOPMENT OF THE PRODUCTION TERM 

The net production term is defined as the difference 
between ionization and recombination rates. Whenever a neu- 
tral atom or molecule loses one or more of its electrons, 
ionization has ocurred. Recombination is the process wher- 
ein positive ions and electrons combine to form neutral 
atoms or molecules. It is simply reverse ionization. If 
the net production is zero then ionization and recombination 
are equal. The net production term is the right side of 
Sgs. 2 and 3, the expression for which is developed here. 

This work considers only three types of recomb ination 

( 1 ) Dissociative - the electron recombines with a molecu- 

lar ion (AB)+, and the recombination energy goes 
into dissociating the molecule and increasing the 
kinetic energy of the resulting products. The pro- 
cess is described as; 

e + (A3) + A ♦ 3 

(2) Three- body (elect ron) - the energy released in the 
recombination process is carried off by a third body 
which is an electron. The process is described as; 

9 ♦ A* e . ' > » A ♦ e 
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(3) Three-body (heav^t^ - the recombina tioa energy is 

carried off by a third body which is a heavy parti- 
cle such as an atom or molecule. The process is 
described as; 

e ♦ + B — > A ♦ B 

This study does not consider other recombination mechan- 
isms such as radiative or dielectric recombination, since 
they are generally of low significance in relation to the 
just described processes. The ions are limited to singly 
positive ions. No account is made for intermediate or 
excited states as this requires a detailed balancing which 
is beyond the scope of this work. Additionally, no accoun-t 
is made of attachment. 

Rewriting Eg. 4 for simplicity of discussion; 

where each refers to the recombination coefficients for 

dissociation, three-body (electron) and three-body (heavy) , 
respectively. 

The recombination coefficients depend upon the ionic 
species involved. some general but representative values 
[Ref. 22 ] follow; 

= 10“ m3/sec 
= 9x1 m3/sec 

= 2x10-13 m3/sec (air) 
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Hinnov and Hirschberg [ Hef. 23] developed a relation for 
three- body (electron) recombination coefficients as a func- 
tion of electron temperature (T^) and density, valid for 
energies less than 0.25 eV. 

1 . 09x10"2On^>r^-9/2 m^/sec 

where the electron temperature is ®K. Gurevich and Pitaen- 
ski [Ref. 24] extended this development to obtain; 

2. 13 x 10 - 2 On^Tg- 9 / 5 in^ m^/sec 

where; 

lnwV= 1.24x102 (T^3/n^)-s 

At 1.0 eV these coefficients are 5 . 6 x 10 * 2 i and 1.05x10~io 
m3/sec respectively. A simplification of the later expres- 
sion for is valid from energies of 0.1 to 10.0 eV and 
n^=10-i« m-3; 

OC-i =0.0 32xT-*- 3 

3rown [ Eef . 18: p.197] gives a value for the dissocia- 
tive recombination coefficient; 

= 1. 5x10-12 m3/s 

The ionization coefficient can be determined from con- 
sidering that there is no net production in the free stream, 
i. e. , far from the electrode. Here; 

0 =>1n -o(,N 2 -A.jN2n^ (B2) 
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where is the stationary state. Solving for and sub- 

stituting into Bg 31 gives; 



n£^(N-n- ) («, + o^no) 


+ o^’jneT^9/2«(N2 


-n^ nc ) 


(B3) 


= 1 .5x10-12 


(from Ref. 


18) 




(/jn^= 1.09x10-20 


(from Ref. 


23) 




<yjn^= 2.0x10-13 


(from Ref. 


22) 





The net production term is the relation, suitably non 
diaensionalized, that appears in the computer subprogram 
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APPENDIX C 



NOMERICAL SOLUTION METHOD 



A. GENERAL DISCUSSION 



There are two basic approaches used in solving a nonli- 
near system of equations, such as possed in this study. A 
given method can usually be classified as one of the ascent 
(or descent) or more commonly refered to as gradient meth- 
ods, or as one of the Newton methods. The method applied in 
this study and by Dolson [Ref. 12] is one of the Newton 
methods. A brief review of both methods will be presented 
to facilitate development of the particular method employed 
in this study. 

Let the nonlinear system of equations be written as; 



(Cl) 



^>7 I * ^2 1 • • • t ) “ ^ 

or in vector form for convenience; 

f ("x) = 0 (C2) 

In the descent method the problem is possed in such a 
way so as to form a function which is zero at any solution 
and positive otherwise. A typical such function would be; 
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(C3) 



0 (xi,X2 , . . .xn) 



n 




2 



so that the problem is reduced to finding a minimum of 0, 
which also happens to be the Talue zero. In solving prob- 
lems by this method, an initial guess to the solution must 
be provided, call it x<i>. Then a direction d is determined 
so as to point to a better solution, This direction 

is often chosen to be the gradient of 0(x) • The operation 
-grad 0(x) is the direction in which 0(x) decreases most, 
and is locally the best direction in which to reduce 0(x) . 

d = -grad 0(x) 



so that ; 



x<2) = z< 1 ) - a*grad 0(x<'>) 

where '‘a'* is chosen to be an appropriate step size. A prob- 
lem is encountered at local minima of 0(x) , which needs to 
be avoided; One must also use the step ”a'* in such a manner 
as to control the acceleration along slowly or rapidly vary- 
ing 0(1). It nay be that by the nature of 0(x) , the gra- 
dient may have to be approximated by finite differences, an 
approach which offers a flexible recourse even when not 
reguired. 
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The Newton method involves approximating the nonlinear 
problem by a suitable linear problem involving derivatives 
of the nonlinear equations evaluated at The linear- 

ized problem is simply a set of linear algebraic equations, 
the solution to which furnishes an increment that when added 
to the solution will yield x<2). This method is a gen- 

eralization to n- dimensions of Newton’s method for evaluat- 
ing the zeros of a function. 

Each function fi (xj is approximated by the initial two 
terms of a Taylor series expansion about the solution vector 



The right-hand-side of equation is the linear vector 
function (tangent hyperplane) of x which best approximates 
the nonlinear function f evaluated at x^ ^ . Se set this 
right-hand-side equal to zero and solve for x, which is an 
improvement from x< ^ . we call this improved solution vec- 
tor ^ , and from equation C4 we get; 



x< I > ; 



fi (^) = fi(x<i>) t J(x<i>)*[x- x<i>] 



(C4) 



where "J " is the familiar Jacobian matrix; 




n 



(C5) 



x(2) = kC 1 ) - J(x< I )) -» »f (X) 



(C6) 
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and often is 



The Jacobian need not be known very accurately, 
left unchanged for several iterations since its evaluation 
often entails penalizing computer time. Thus the method 
(slightly modified as described briefly in section III.B) 
used by Dolson [Ref. 12], is referred to in many texts as 
the Newton method or as the Newton-Raphson method [Refs. 25, 
26 ]. 

The Jacobian matrix of coefficients for elliptic prob- 
lems is typically a sparse banded matrix, to which there are 
many computer subroutines available for solution. 

B. MODIFIED NEHTON-RAPHSON METHOD 

The particular solution technique employed in this study 
merits a detailed discussion, since it is not common in the 
literature. Bailey and Touryan [Ref. 27] use a similar 
method and provide some description of its development. The 
method is similar to other Newton methods however, the Jaco- 
bian matrix is replaced by a matrix of coefficients devel- 
oped from a sequence of linear differential equations which 
approximate the system of nonlinear equations. This method 
will be referred to as a modified Newton-Raphson method or 
MNH method for short, and will be presented in general using 
a specific example. 
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Take Poisson' s equation given earlier and repeated here 
for simplicity; 

V20= =, (oe -°t) (C7) 

The known solution variables are; 

and n^ < » > 

In general the first st^ is linearization of the equa- 
tion. In this case Sq. C7 is linear, but it will serve as 
an example. ie assume that the suceeding solution variables 
at iteration- (2) , namely z<2>, can be approximated by the 
solution zCi) plus some small change ^z<‘>. 

^2< 2> = 2 < 1) ♦ Z< 1 ) (C8) 

where; 

z<k) = 

It remains to determine 4z^i^. Substitution of this 
relationship into the Poisson Equation is the second step ; 

C 7VV C,(n€”- -V 

rearranged ; 

7^V'-C,(At'e’ - AHi'”) = - [vV'- C, ( He“'- OJ (C9) 
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The right hand side of Eg. C9 is a function of the known 
solution vector z<i>. We are applying this relationship at 
each computational node, k second order partial differen- 
tial eguation can be approximated by a second order differ- 
ence equation or ’*3-point formula," so Eg. C9 will become a 
system of linear algebraic relationships involving near 
neighboring values of the solution vector z<‘>. The left 
side of Eg. C9 is composed of a system of coefficients 
expressible in terns of z<i) times the unknown (perturba- 
tion) vector ^z<»>. So we now have; 

A(z< »))• »> = -F(z<‘>) (CIO) 

Here the matrix of coefficients "A" is a sparse matrix 
with a relatively narrow band of non-zero elements about the 
main diagonal, provided that the "natural" order of the 
resulting linearized equations is preserved. There are 
standard computer routines available which take advantage of 
this banded structure, such as LEQTIB, an IMSL routine. 

It follows that since ^z<»> = z^*) - z<‘>; 

Z<2) = Z<1> - C A(Z< 1>) (Z<‘ >) (C11) 
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For the case of a nonlinear term such as; 



there is no additional difficulty. When linearized this 
term becomes; 

^ }fU 4. (Cl 2) 

Jx <5x ^ 

We neglect high order terms created by the inner products of 
the elements of such as the last term in relation 

C12. The first term would be a part of the function F(z<i>) 
in the right side of Equation CIO , and the remaining two 
terms would be included as part of the expression: 

A (Z< ‘ >) • ^ Z< i >. 
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APPENDIX D 



DESCRIPTION OP COMPUTER PROGRAMS 

This appendix serves as a user oriented guide to the 
computer prograas listed in App. E. The programs are listed 
in the order used to generate solutions to the problem for- 
mulated in this report. The output consists of printed 
results for each iteration and punched output to provide a 
permanent file copy, which is to be used as input to the 
three plotting programs also listed in App. E. These plot- 
ting programs use the Vecsatec plotter and associated soft- 
ware, and represent the only facility not commonly available 
elsewhere. The main program SHTH as presented, is config- 
ured to also produce output at a monitor terminal as used on 
CP/cas (laa 360) or avs (IBH 370). The output in summary is 
via three devices; 

DEVICE 4 - typically a CRT 

DEVICE 6 - typically a printer 

DEVICE 7 - typically a punch or disk output 

These can be manipulated in convenient configurations to 

give a full or partial (sampled) display of the solution at 

each iteration, or merely to display convergence information 

in a simple format. - 



■) 
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The SHTH program is the main program which together with 
its 9 subroutines and 1 IMSL subroutine solves the non-di- 
mensionalized forms of Eqs. 1 through 4. The input to the 
program consists of control cards and data cards, namely; 

1 control card 

3 or 104 initial solution data cards 

1 control card 



The control card is described here in terms of the Fortran 



variable names thereon. 



ITHIH - IS the iteration counter (ITBH) starting value. 

If ITMIN = 1, then the program will read in the 
next 3 cards as data which contain the short form 
starter solutions. If ITHIH is greater than 1, 
then the program will read in the next 104 cards as 
data, representing a detailed starting solution. 
This detailed starting solution is typically the 
output of a previous run of the program. If ITHIH 
= 0.0, then the program will not read in any data, 
being thus informed that the problem is over. 

ITHAX - is the last iteration value. The oroaram ceases 
iterations at this point and outputs' the solutions 
and auxiliary information. 

ITHDT - enables the employment of the net production term 
(RTSIDE) when ITEB is greater than ITHDT, otherwise 
the net production is set to zero. 

ITEflP - for IT?R less than ITEMP, the electron tempera- 
ture is isothermal, otherwise it is a function of 
S/n. (see SHTH statements 899-990) . 

ITBC - for ITER greater than 5 and less than ITBC, the 
wall boundary conditions for the species is 
floated. 

LADJ - 1 electron temperature is calculated 
- 0 electron temperature is isothermal. 

LPOl - 1 output is punched on device 7. - 0 output is 

not punched. 

H - sets the mesh spacing egual to H'^LAMS (sheath 
length) . 

PKTR - is the adjustment factor to the temperature equa- 
tion (see Sec. 17) and is used when LADJ =1 and . 
ITER IS less than ITEHP. The adjustment factor is 
not used when LADJ is set to zero. When the temp- 
erature is coupled and LADJ =1, then the program 
calculates its own adjustment factor. 

HAME - is reserved for user labelling and identification 
which is used at output. 
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The solution set is; 

PH voltage potential (volts) . 

HE electron density (10'* m“*) 

HI ion density (10'® m”*) 

TH nondioens ional electron temperature 
The short form starter solution is illustrated in App. 

E. It consists of PH, HE and HI values at the anode, wall 
and equilibrium boundaries. The long form starting solution 
consists of PH, HE, HI, and TH distributions, each one 
represented by 25 data cards and a label card. Each card 
has 7 values, (see format-803). 

A typical iteration series would sta^t out for about ten 
iterations from the short form data with ten volts at the 
anode and other appropriate boundary conditions. The ini- 
tial solution converges very well for H = 0.5, the net pro- 
duction term coupled after 2 or 3 iterations, and the wall 
conditions not floated. This output is then used as input. 
The voltage is amplified to say 25 or 35 volts and a new 
control card is structured to couple the temperature or wha- 
tever else is desired. The H value may also be changed to 
reflect a smaller mesh spacing. It must be Jcept in mind 
that too many changes in one solution run can be catas- 
trophic. 
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Each iteration of the solution by the SflTH program 
accomplishes the following; 

(1) floats the NE, NI wall values 

(2) sweeps through the A-matrix and C-vector setting 
up the equations that satisfy the boundary condi- 
tions and the internal mesh conditions 

(3) couples the electron temperature equation 

(4) couples the net production term 

(5) calculates the anode electron density required to 
satisfy current matching btween the anode and the 
f ree-stream. 

The system of equations has the form presented in Eq. 

CIO in App. C, and is solved by use of the IMSL routine 
LEQT1B. The solution vector is then updated from the return 
information contained in vector C. 

There are two solution convergence monitors. The first 
is the sun of the squares of the differences (spectral norm) 
of the -F(z<^>) values from Eq. CIO. The other is the spec- 
tral norm of the upon return from LEQT1B. Both of 

these are available at the video terminal (device 4) and 
printer (device 6,) . 

After solution update, the elctron temperature is calcu- 
lated and adjusted per the control card options. It should 
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be pointed out that on the last iteration, the temperature 
distribution is automatically calculated for output, but is 
not used in the scheme, unless otherwise coupled per the 
control card. 

The matching of the anode current to the free-stream 
current is accomplished by solving for the anode electron 
density using subroutine INODE. The call to this subprogram 
will require a change to the SHTH program as it is presently 
overridden with a comment card, (prior to statement 9000) . 

The final outputs are straight-forward. Subroutines 
J8ALL and JOOLE are called in this phase of the program, in 
order to compute the current into the wail and the Joule 
heat distribution. The electric field, net production and 
space charge density distributions are also calculated and 
provided as punched and printed output. 

The punched output is used as input for the STREAM, PICT 
and PLOT programs. These programs use the 7ersatec plotter 
to provide a presentation of the current stream lines, an 
oblique perspectives of the solutions and plots of the solu- 
tions as they vary normal to the wall extending from the 
anode site. The details of these programs are combersome, 
and it is wise to leave them to the user for inspection. 
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The first two distributions in the punched output are 
used in the STBEil! program, requiring that the nonpertinent • 
data be removed. The next 8 distributions are used for the 
PLOT and PICT programs. The PICT program has an option card 
filled with ones or zeros in the first eight positions. 

This option card enables the user to elect which distribu- 
tions to plot. The plots provided in this report as Figs. 

12 through 14, were produced using the PICT and PLOT pro- 
grams. The programs are simple and speak for themselves. 
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APPENDIX E 



COMPOTEE PROGRAM LISTINGS 

This appendix contains a listing of the programs used in 
the solution of the equations presented in this work 
(namely: Bqs. 9,10 and 11), and the auxiliary programs that 
provided the plots for the report. 

SHTH (MAIN PROGRAM) 109 

DELOP 127 

GEDNT 129 

TEMP 131 

3TSIDE 132 

PATCH 133 

ANODE * 134 

JHALL 136 

JOULE 137 

OUTPUT 138 

PLOT 139 

PICT 141 

STREAM 144 

INITIAL INPUT (sample) 150 
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